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Abstract.  An  important  problem  for  space  vehicles  is  contamination  of  the  spacecraft  surface  by  combustion  products  and 
unburned  fuel  exhausting  from  control  thruster  nozzles.  The  objective  of  the  present  work  is  a  detailed  simulation  of  the 
backflow  in  the  case  of  gas  exhaustion  from  the  nozzle  into  vacuum.  A  specific  feature  of  this  problem  is  the  drastic  expansion 
of  the  initially  dense  flow  at  the  nozzle  exit,  which  makes  necessary  simulations  of  the  near-continuum  flow  inside  the  nozzle, 
the  transitional  flow  near  the  nozzle  exit,  and  the  free-molecular  flow  outside  the  nozzle.  The  simulations  are  performed  by  a 
combined  Navier-Stokes/DSMC  approach. 


INTRODUCTION 

An  important  problem  arising  in  long-term  operation  of  space  vehicles  is  contamination  of  the  spacecraft  surfaces 
with  the  combustion  products  issuing  from  the  nozzles  of  the  control  thrusters  (see  ref.  [1]).  During  thruster  firing,  a 
large  amount  of  propellant  is  exhausted  into  space,  and  droplets  forming  as  a  result  of  liquid  fuel  condensation  and 
breakup  of  the  nozzle  cooling  film  can  reach  the  spacecraft  surface.  One  of  the  possible  mechanisms  provoking  the 
reverse  motion  of  liquid  droplets  may  be  the  carrier-gas  backflow  originating  from  the  nozzle  boundary  layer,  which 
turns  sharply  around  the  nozzle  lip  while  undergoing  continuum  expansion  and  then  passing  into  transitional  and  free- 
molecular  regimes  as  the  density  decreases.  It  can  be  assumed  that  the  gas  flow  drags  the  droplets  from  the  nozzle  into 
the  backflow  region. 

In  this  paper  we  simulate  the  gas  flow  exhausting  from  a  supersonic  nozzle  with  Reynolds  number  Re  =  590.  In 
addition,  we  performed  a  study  of  the  mechanisms  of  the  flow  turn  around  the  nozzle  lip  and  the  backflow  origination. 
A  specific  feature  of  this  problem  is  the  drastic  expansion  of  the  initially  dense  flow,  which  makes  necessary  using  the 
combined  Navier-Stokes/DSMC  approach. 


NUMERICAL  APPROACH 

The  two-dimensional  and  axisymmetric  Navier-Stokes  equations  are  solved  with  a  shock-capturing  MUSCL  TVD 
scheme.  The  standard  van  Leer’s  formula  with  the  minmod  slope  limiter  [2]  is  used  for  reconstruction  of  cell-faced 
variables  from  cell-centered  ones.  The  viscous  terms  are  approximated  by  the  central  finite  differences.  The  HLLEM 
approximate  Riemann  solver  [3]  is  employed  to  calculate  the  inviscid  fluxes.  The  HLLEM  solver  is  constructed  to  be 
robust  near  low  densities  and,  therefore,  it  is  particularly  appropriate  for  simulation  of  rapidly  expanding  flows  such 
as  in  the  present  paper. 

For  low-Reynolds-number  regimes,  the  initial  rarefaction  effects  are  taken  into  account  by  means  of  imposing 
velocity  slip  and  temperature  jump  boundary  conditions  on  solid  walls.  For  this  purpose,  the  slip  conditions  deduced 
by  Kogan  [4]  are  used. 

The  DSMC  code  SMILE  [5]  developed  in  Computational  Aerodynamics  Laboratory,  Institute  of  Theoretical  and 
Applied  Mechanics  is  utilized  for  kinetic  modeling.  Intermolecular  collisions  are  computed  by  the  Variable  Hard 
Sphere  (VHS)  model.  The  number  of  collisions  is  calculated  by  the  majorant  frequency  method.  In  the  case  of 
a  polyatomic  gas,  the  energy  redistribution  between  rotational  and  translational  degrees  of  freedom  is  obtained  in 
accordance  with  the  Larsen-Borgnakke  model  [6]  with  a  constant  rotational  collision  number.  Up  to  9.6T06  molecules 
are  used  to  simulate  the  flows  under  consideration.  The  computations  are  performed  on  a  32-processor  Linux  cluster. 
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FIGURE  1. 


An  example  of  the  hybrid  continuum/kinetic  simulation  of  low-Reynolds-number  flow  in  a  conical  nozzle. 


In  the  course  of  hybrid  simulations,  a  coupling  of  Navier-Stokes  solver  and  the  DSMC  method  is  based  on  the 
use  of  the  NS  solver  in  the  nozzle  and  in  the  plume  core  and  the  DSMC  method  in  the  neighborhood  of  the  nozzle 
lip  and  in  the  backflow  region.  Flow  parameters  on  the  boundaries  of  the  DSMC  computational  domain,  which  are 
located  inside  the  nozzle  and  in  the  central  part  of  nozzle  plume,  are  taken  from  a  previous  Navier-Stokes  simulation. 
The  velocity  distribution  function,  which  is  needed  to  provide  necessary  boundary  conditions  for  the  DSMC  method, 
is  specified  as  the  local  Maxwellian  distribution.  It  is  worth  noting  that,  from  the  general  viewpoint,  the  Chapman- 
Enskog  distribution  function  is  a  natural  choice  for  model  particles  coming  into  the  DSMC  computational  domain 
at  the  boundary  between  the  NS  and  DSMC  domains.  However,  the  Chapman-Enskog  distribution  function  has  a 
drawback:  it  becomes  negative  with  deviation  from  the  Maxwellian  distribution.  Moreover,  it  was  demonstrated  in  a 
number  of  works  (see,  e.g.  [7])  that  the  Chapman-Enskog  distribution  function  may  not  provide  much  improvements 
against  the  Maxwellian  distribution. 


RESULTS 

Gas  exhausting  from  nozzle  was  modelled  on  the  first  stage  of  computations.  The  nozzle  geometry  and  flow  parameters 
were  taken  from  Rothe’s  experiments  [8].  Corresponding  nozzle  and  flow  parameters  are:  throat  radius  K,  =  2.55  mm, 
exit  radius  Re  =  20.92  mm,  stagnation  temperature  7b  =  300  K,  stagnation  pressure  po  =  945  Pa,  Reynolds  number 
Re  =  590.  The  density  flowfields  computed  in  this  case  with  the  Navier-Stokes  solver  inside  the  nozzle  and  by  the 
DSMC  method  for  the  external  flow  are  shown  in  Fig.  1.  The  domains  of  the  computations  made  by  these  two  methods 
overlap  in  some  region  near  the  nozzle  exit.  The  boundaries  of  computational  domains  are  shown  with  red  lines. 
The  density  contours  are  shown  as  solid  (for  Navier-Stokes)  and  dashed  (for  DSMC)  curves.  From  this  figure,  close 
agreement  between  two  solutions  in  the  common  part  of  two  computational  domains  is  evident. 

The  density  distribution  along  the  nozzle  axis  and  the  density  distribution  across  the  nozzle  close  to  the  nozzle 
exit  are  shown  in  Fig.  2.  Here  the  results  of  the  Navier-Stokes  simulations  (with  slip  and  no-slip  conditions  on  the 
wall)  and  the  DSMC  results  are  compared  with  the  experimental  data  obtained  in  Rothe’s  experiments.  The  numerical 
results  along  the  axis  compare  well  with  the  experiment  along  the  whole  nozzle  length.  However,  as  is  seen  in  Fig.  2, 
the  flow  near  the  nozzle  axis  is  poorly  reproduced  in  Navier-Stokes  computations  with  the  no-slip  conditions.  The  slip 
conditions  yield  results  which  are  closer  to  the  experiment. 

In  addition,  we  investigated  the  flow  around  the  nozzle  lip  for  three  nozzle  lip  shapes:  sharp  lip,  and  two  rounded 
lips  with  the  radii  1.67  •  1 0  4  m  and  6. 1 1  •  1 0  4  m.  The  DSMC  simulations  were  performed  in  a  small  domain  near  the 
nozzle  lip.  The  inflow  data  were  input  into  the  DSMC  simulations  from  the  results  of  the  Navier-Stokes  computations. 

Figure  3  shows  the  mass  flux  isolines  for  three  different  nozzle  lip  shapes.  Integral  mass  flow  rate  in  the  reverse 
direction  normalized  on  total  mass  flow  rate  through  the  nozzle  was  1.46-  10~3  for  the  nozzle  with  a  sharp  lip, 
1.33  •  1 0  3  -  for  the  nozzle  lip  rounded  with  the  radius  1.67  •  1 0  4  m,  and  9.91  •  1 0  4  -  for  the  nozzle  lip  rounded 
with  the  radius  6.11-  1 0  4  m.  It  is  therefore  evident  that  the  higher  is  the  rounding  radius  of  the  nozzle  lip  the  lower  is 
the  backflow  component  of  the  mass  flow  rate  from  the  nozzle. 

To  further  investigate  the  origin  of  the  backflow  we  performed  additional  computations  of  the  flow  exhausting  from 
the  nozzle.  Outside  the  nozzle,  in  the  backflow  region,  we  set  up  a  screen  located  at  some  height  above  the  outer 
contour  of  the  nozzle  (see  Fig.  4).  The  purpose  of  this  screen  was  to  separate  different  backflow  components:  let  the 
backflow  originated  near  the  nozzle  wall  go  through,  but  eliminate  possible  backflow  from  the  periphery  of  the  nozzle 
plume. 
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FIGURE  2.  Density  distribution  along  the  nozzle  axis  (left)  and  across  the  nozzle  at  X /R*  =  18.7  (right) 


FIGURE  3.  Mass  flux  fields  for  three  different  shapes  of  the  nozzle  lip 


Figure  4  presents  density  isolines  in  the  backflow  region  obtained  with  and  without  the  screen.  As  seen  in  the  figure, 
the  density  flowfield  is  changed  at  some  distance  from  the  outer  contour  of  the  nozzle  and  is  unchanged  near  the  outer 
nozzle  wall.  The  density  and  mass  flux  distributions  in  radial  direction  behind  the  screen  are  shown  in  Fig  5.  Is  is  clear 
that  the  screen  The  presence  of  the  screen  leads  to  a  very  significant  decrease  in  density  on  the  left  side  of  the  screen.  It 
means  that  the  most  part  of  the  backflow  gas  comes  from  the  plume  rather  than  from  the  nozzle  boundary  layer.  It  can 
be  assumed  that  the  mechanism  of  backflow  origination  from  the  plume  is  purely  kinetic  —  molecules  on  the  plume 
periphery  gain  an  oppositely  directed  velocity  owing  to  intermolecular  collisions. 


FIGURE  4.  Backflow  formation  on  the  lip:  without  screen  (left)  and  with  screen  (right) 
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FIGURE  5.  Radial  distribution  of  density  (left)  and  mass  flux  (right)  behind  the  screen 

VISCOUS  FLOW  AROUND  AN  EXPANSION  CORNER 

Based  on  the  results  described  above,  it  is  evident  that  the  investigation  of  flow  expansion  around  the  nozzle  lip  is 
critically  important  for  studying  backflow  contamination.  For  an  inviscid  flow,  it  is  well  known  from  the  Prandtl-Meyer 
solution  that  there  is  a  limiting  flow  deflection  angle.  At  this  limiting  angle,  all  internal  energy  of  the  gas  transforms 
into  kinetic  energy  of  motion,  so  that  the  flow  temperature  becomes  zero  while  the  Mach  number  keeps  growing  up 
to  infinity.  The  viscosity  effects  can  change  the  expansion  flow  substantially  and,  in  particular,  they  favor  flow  turning 
by  large  angles  and  backflow  formation.  To  investigate  this  phenomenon  in  more  detail,  we  consider  here  a  model 
problem  consisting  of  the  two-dimensional  flow  over  an  expansion  corner  formed  by  a  plane  wall  (flat  plate)  suddenly 
deflecting  by  a  large  angle.  The  wall  deflection  angle  (pw  is  varied  in  different  simulations.  To  exclude  any  influence 
of  models  governing  the  exchange  between  translational  and  internal  degrees  of  freedom  of  molecules,  a  monatomic 
gas  (y  =  5/3,  Argon)  is  considered.  The  simulations  were  performed  with  inviscid  conditions  (non-permeability  for 
Navier-Stokes  simulations,  specular  reflection  for  DSMC).  The  simulations  were  performed  for  the  following  flow 
parameters:  the  inflow  Mach  number  =  4.38  and  static  flow  temperature  7’00=75  K.  The  unit  Reynolds  number 
was  varied  Re  =  476,4760,47600  to  investigate  the  effects  of  flow  viscosity.  The  flow  was  simulated  by  the  DSMC 
method  and  with  the  Navier-Stokes  equations.  The  Navier-Stokes  computations  were  performed  for  the  expansion 
corner  of  (pH.  =  37.1°,  which  coincides  with  the  maximum  Prandtl-Meyer  angle  <pn]ax  for  y  =  5/3  and  M„  =  4.38. 
Note,  the  maximum  Prandtl-Meyer  angle  for  y  =  5/3  and  M„o  =  1  is  90°.  The  DSMC  computations  were  performed 
with  (pn,  =  37.1°  and  also  with  <pM,  =  180°,  the  latter  simulates  expansion  into  vacuum  around  the  nozzle  lip.  The  flow 
was  assumed  to  be  two-dimensional. 

The  typical  flowfields  obtained  are  shown  in  Fig.  6.  These  computations  were  performed  at  the  Reynolds  number 
Re  =  47600,  the  wall  defelection  angle  in  the  DSMC  simulation  was  <pH  =  180°.  A  close  resemblance  of  two  flowfields 
is  worth  noting,  it  holds  in  spite  of  the  different  shape  of  the  computational  domains  used.  Further,  it  is  evident  from 
the  flowfields  that  the  numerical  solutions  are  close  to  the  classical  Prandtl-Meyer  solution  only  in  some  range  of  polar 
angles  <p  (here  (p  =  0  corresponds  to  the  axis  going  along  Wall  1,  and  (p  is  positive  in  the  clockwise  direction).  Another 
interesting  feature  of  the  solutions  obtained  is  a  non-monotonous  behavior  of  the  Mach  number:  starting  from  some  <p 
it  begins  to  increase. 

The  angular  distributions  of  the  flow  parameters  (non-dimensional  density  p/p0 «,  Mach  number,  temperature  T /7C 
and  the  flow  velocity  normalized  on  maximum  adiabatic  velocity  \ii\/um)  are  given  in  Figs.  7  and  8.  These  distributions 
were  recorded  along  the  circle  with  radius  r/L  =  0.4  centered  at  the  vertex  of  the  expansion  corner.  For  comparison, 
the  inviscid  Prandtl-Meyer  analytical  solution  is  also  plotted  in  the  figures.  The  DSMC  computations  in  this  case  were 
performed  for  wall  deflection  <p„.  =  <pmax  =  37.1°  and  also  for  flow  expansion  into  vacuum  around  the  trailing  edge, 
i.e.  for  wall  deflection  <p„.  =  180°.  These  results  are  very  close  to  each  other,  hence  only  (pw  =  <pmax  =  37.1°  data  are 
given  in  the  figures. 

It  is  evident  that  the  flow  viscosity  substantially  changes  the  flow  around  the  corner  even  if  the  inviscid  wall  boundary 
conditions  are  imposed.  This  is  manifested  in  termination  of  growth  of  the  flow  Mach  number  at  some  polar  angle, 
whose  value  depends  on  the  Reynolds  number,  and  also  in  an  excessive  increase  in  flow  temperature  compared  to  the 
Prandtl-Meyer  solution.  These  changes  become  more  apparent  as  the  flow  Reynolds  number  decreases.  Note,  the  flow 
parameters  tend  to  their  inviscid  limit  with  increasing  Reynolds  number.  It  should  be  noted  that  <pmax  is  based  on  the 
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FIGURE  6.  Mach  flowfields  computed  by  the  DSMC  method  (upper)  and  by  the  Navier-Stokes  solver  (lower) 


FIGURE  7.  Angular  distribution  of  density  (left)  and  Mach  number  (right)  around  the  expansion  corner 


FIGURE  8.  Angular  distribution  of  temperature  (left)  and  velocity  (right)  around  the  expansion  comer 
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inviscid  theory  and  does  not  necessarily  hold  in  this  case  where  the  viscous  forces  produce  their  effects.  Still,  in  the 
absence  of  the  boundary  layers  on  the  walls  that  form  the  expansion  corner,  the  inviscid  theory  predicts  the  maximum 
turning  angle  quite  well.  Nevertheless,  the  density  in  the  Navier-Stokes  computations  performed  with  Re  =  476  is 
quite  different  from  zero  even  for  <p  =  <pmax .  The  latter.  Re  =  476,  case  demonstrates  the  most  prominent  differences 
of  the  Navier-Stokes  results  compared  to  the  DSMC  computations.  These  differences  are  evident  in  this  case  from  the 
very  beginning  of  flow  expansion.  The  flow  Mach  number  in  Navier-Stokes  computations  does  not  grow  at  all.  This 
is  caused  by  the  growth  of  flow  temperature  due  to  the  work  of  viscous  forces.  Based  on  this,  we  can  conclude  that 
the  Navier-Stokes  equations  are  not  applicable  to  simulation  of  flow  expansion  at  such  a  low  Reynolds  number.  For 
higher  Reynolds  numbers,  i.e.,  for  Re  =  4760  and  Re  =  47600,  the  Navier-Stokes  results  are  in  close  agreement  with 
the  DSMC  computations  up  to  the  polar  angles  <p  «  5°  and  <p  ~  15°,  respectively. 

We  can  conclude  that  the  use  of  the  combined  Navier— Stokes/DSMC  approach  for  simulation  of  the  flow  around 
the  nozzle  lip  is  essentially  important  since  the  Navier-Stokes  approach  become  no  longer  acceptable  as  the  flow  is 
expanded. 


CONCLUSIONS 

A  detailed  numerical  investigation  of  the  backflow  formation  phenomenon  for  nozzle  plumes  expanding  into  vacuum 
has  been  performed  by  a  hybrid  approach  combining  the  continuum  (Navier-Stokes)  and  kinetic  (DSMC)  computa¬ 
tions. 

The  flow  around  the  nozzle  lip  has  been  investigated  in  detail  for  different  shapes  of  the  nozzle  lip  (sharp  and 
rounded).  It  has  been  found  that  the  mass  flow  rate  in  the  backflow  region  is  rather  small,  amounting  approximately  to 
1 0  3  of  the  total  mass  flux  through  the  nozzle.  The  simulation  also  demonstrate  that  intermolecular  collisions  in  the 
plume  peripheral  zone  contribute  substantially  to  formation  of  the  low-density  backflow. 

The  viscous  expansion  flow  over  an  expansion  corner  has  been  investigated  to  elucidate  details  of  backflow 
formation  near  the  nozzle  lip.  The  results  of  Navier— Stokes  and  DSMC  simulation  have  been  compared  and  analyzed. 
It  has  been  found  that  the  effects  of  viscosity  become  more  significant  with  the  flow  expansion  and  dominate  in 
the  rarefied  regime.  The  influence  of  viscous  dissipation  causes  a  decrease  in  the  Mach  number  and  a  growth  of 
temperature  in  contrast  with  the  inviscid  Prandtl-Meyer  solution. 
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